function    A210()
format long;

% Monte Carlo法，一种基于随机数的数值法，通过对各影响量进行大量离散抽样进行数值计算来解决不确定度评定的问题。
% 测量模型的函数关系
tic                                       % 计时开始
% 参数设置
num_points = 1000; % 数据点数量
rayleigh_scale = 0.5; % 瑞丽分布的尺度参数
normal_mean = 0; % 正态分布的均值
normal_std = 0.2; % 正态分布的标准差
mix_ratio = 0.7; % 瑞丽分布与正态分布的混合比例（0~1之间）
% 生成瑞丽分布数据
rayleigh_data = raylrnd(rayleigh_scale, num_points, 1);
% 生成正态分布数据
normal_data = normrnd(normal_mean, normal_std, num_points, 1);
% 数据混合
out_y = mix_ratio * rayleigh_data + (1 - mix_ratio) * normal_data;
% 输出量的可视化
subplot(2, 2, 1)
set(gca, 'Fontsize', 12);
plot(1:length(out_y), out_y, 'k-');
title('输出量分布');
hold off

[PDF_y, CDF_y] = ksdensity(out_y);        % 计算输出量概率密度估计值
%PDF_y_normalized = PDF_y / max(PDF_y);
%CDF_y_normalized = CDF_y / max(CDF_y);
% 计算测量数据的矩(mean，std，skewness，kurtosis)
mu1 = mean(out_y);
sigma1 = std(out_y);
skew1 = skewness(out_y);
kurt1 = kurtosis(out_y);
% 输出量概率密度曲线的可视化
subplot(2, 2, 2)
set(gca, 'Fontsize', 12);
hold on
%plot(1:length(PDF_y), PDF_y_normalized, 'r-');
%plot(1:length(CDF_y), CDF_y_normalized, 'b-');
plot(CDF_y, PDF_y, 'r-');
xlabel('输出量');
ylabel('概率密度');
title('输出量概率密度分布');
legend(['均值（', num2str(mu1), '），方差（', num2str(sigma1), '），偏态（', num2str(skew1), '），峰态（', num2str(kurt1), '）']);
hold off

% 包含区间
sort_y2 = sort(out_y);
for m = 1 : 1000
    Z2(:, m) = prctile(sort_y2, m/10);
end
for k4 = 951 : 1000
    in3_X2(:, k4-950) = Z2(k4)-Z2(k4-950);
end
J2 = find(in3_X2==min(in3_X2));                                  % 计算最短包含区间所对应的位置
Y2 = [mean(out_y), std(out_y), Z2(J2), Z2(951+J2)];      % 给出平均值、标准差、95%概率包含区间的左端点和右端点
% 区间可视化
subplot(2, 2, 3)
histogram(out_y, 100)
hold on
a2 = 0: 1 : 35;
for i = 1 : length(a2)
    per_mu2(i) = Y2(1);
    per_left2(i) = Y2(3);
    per_right2(i) = Y2(4);
end
plot(per_mu2, a2, 'r-')
plot(per_left2, a2, 'r--')
plot(per_right2, a2, 'r:')
xlabel('输出量');
ylabel('概率密度');
title('输出量的区间');
legend('数据直方图', ['均值', num2str(Y2(1))], ['左端点', num2str(Y2(3))], ['右端点', num2str(Y2(4))])
hold off

% 使用 out_y 数据作为蒙特卡洛输入
num_samples = length(out_y); % 保持样本数量一致
J4 = 1000; % 蒙特卡洛模拟次数

% 初始化蒙特卡洛模拟结果矩阵
model_X4 = zeros(J4, num_samples);

for h4 = 1:J4
    % 在 out_y 基础上引入随机扰动
    perturbation = normrnd(0, 0.05, size(out_y)); % 添加小扰动
    model_X4(h4, :) = out_y(:)' + perturbation'; % 确保大小一致
end

% 计算蒙特卡洛模拟的统计结果
mu_Y4 = mean(model_X4, 2); % 每次模拟的均值
std_Y4 = std(model_X4, 0, 2); % 每次模拟的标准差
left_Y4 = prctile(model_X4, 2.5, 2); % 95% 包含区间左端点
right_Y4 = prctile(model_X4, 97.5, 2); % 95% 包含区间右端点

% 聚合所有模拟结果
model_X4_all = model_X4(:);
mu_out = mean(model_X4_all);
std_out = std(model_X4_all);
left_out = prctile(model_X4_all, 2.5);
right_out = prctile(model_X4_all, 97.5);

% 偏差计算
dlow = abs(mu_out - left_out);
dhigh = abs(mu_out - right_out);
subplot(2, 2, 4);
hold on;

% 获取直方图数据（原始计数）
hist_counts = histogram(out_y, 100);
bin_counts = hist_counts.Values;   % 每个柱的计数
bin_edges = hist_counts.BinEdges; % 边界
bin_centers = (bin_edges(1:end-1) + bin_edges(2:end)) / 2; % 中心位置
delete(hist_counts); % 删除原始直方图（重新绘制归一化版本）

% 归一化直方图 Y 值
bin_counts_normalized = bin_counts / max(bin_counts); % 最大值归一化
bar(bin_centers, bin_counts_normalized); 

% 绘制 out_y 的概率密度曲线
[f_out, xi_out] = ksdensity(out_y);
f_out = f_out / max(f_out); % 归一化到 [0, 1]
plot(xi_out, f_out, 'b-');

% 计算蒙特卡洛模拟的概率密度，并归一化
[f_MC, xi_MC] = ksdensity(model_X4_all);
f_MC = f_MC / max(f_MC); % 归一化到 [0, 1]

% 绘制蒙特卡洛模拟的概率密度曲线
plot(xi_MC, f_MC, 'r--');

% 绘制蒙特卡洛模拟的包含区间
plot([left_out, left_out], [0, 1], 'r--');
plot([right_out, right_out], [0, 1], 'r--');

% 设置图例和标题
xlabel('输出量');
ylabel('归一化概率密度');
title('MCM与输出数据的概率分布对比');
legend('输出归一化直方图', '输出概率密度', 'MCM概率密度', ...
       'MCM左端点', 'MCM右端点', 'Location', 'Best');
hold off;


toc                                       % 计时结束

